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The author is developing a numerical code with thousands of emission zones to simulate the time-dependent 
multi- waveband emission from blazars. The code is based on a model in which turbulent plasma flowing at a 
relativistic speed down a jet crosses a standing conical coUimation shock that accelerates electrons to maximum 
energies in the 5-100 GeV range. This paper reports early results produced by the model. The simulated light 
curves and time profiles of the degree and position angle of polarization have a number of features in common 
with the observational data of blazars. Maps of the polarized intensity structure can be compared with those 
of blazars observed with very long baseline interferometry at short millimeter wavelengths. 



1. INTRODUCTION 

Reproduction of the time-dependent multi- 
waveband emission of blazars poses severe theoretical 
challenges. Observations with the Ver y Long Baseline 
Array (VLBA) at 7 mm wavelength [Marscher et alJ 
I2OI2] have found that the majority of 7-ray flares in 
blazars are simultaneous with either the passage of a 
new superluminal knot through the bright stationary 
feature termed the "core" or a brightening of the core 
with no prominent knot subsequently appearing. The 
core lies parsecs fro m the central engine in blazars 
[see iMarscheiJ I2OO6I for a discussion]. Nevertheless, 
time-scales of flux changes can be as short as hours, 
even for very high-energy 7 rays that must also be 
emitted on parsec sc ales to avoid a hig h opacity to 
pair-production [e.g. jAleksic et aLl l201l!|. 

One way to explain the rapid variability so far from 
the central engine is to imagine that the emission 
arises from localized regions where Doppler factors 
exceed the values derived from apparent superlumi- 
nal speeds of radio knots. Two possibilities to real- 
ize this s cenario have been pro posed: magnetic recon- 
necti ons iGiannios et al. 20091 and relativistic turbu- 



lence [Marscher and Jorstadl2010l . lNaravan and PiranI 
|2012| |. The author is developing a numerical code, 
called the Turbulent Extreme Multi-Zone (TEMZ) 
model, that can compute the time-dependent spec- 
tral energy distribution (SED) and therefore multi- 
frequency light curve for these physical scenarios. 
Thus far, the code has been applied to an approxi- 
mate simulation of turbulent plasma crossing a stand- 
ing conical shock in the jet. It includes two sources 
of seed photons for inverse Compton scattering: a 
hot, dusty torus and a Mach disk (also often called 
a "working surface"), which is a strong shock ori- 
ented transverse to the jet flow that occupies a small 
portion of the jet cross-sectio n centered on the axis 
Courant and Friedrichsl Il948| | . One of the primary 



merits of this effort is that the simulations can be 
compared with the wealth of data produced by multi- 
waveband monitoring programs: light curves, SEDs, 



and polarization vs. time, frequency, and position. 
This progress report presents early results from this 
effort. 



2. DESCRIPTION OF THE TEMZ MODEL 

The code is based on a model in which turbu- 
lent plasma flowing at a relativistic speed down a jet 
crosses a standing conical coUimation shock that ac- 
celerates electrons to energies up to some maximum 
value, probably in the range from 5 to 100 GeV. The 
code divides the plasma into many cells (see Fig. [T|), 
each with a different direction of magnetic field. This 
provides a crude simulation of turbulence. The en- 
ergy density of electrons is modulated randomly such 
that the variations in flux produce a power spectral 
density described by a power law with a slo pe similar 
to th e observed value in the range of 1.2-2 [Larsson I 
I2OI2J I . The direction of the magnetic field is different 
for each cell, and the maximum electron energy gener- 
ated at the shock can either depend on that direction 
or be set randomly within a probability distribution. 
The plasma in each cell has a turbulent velocity with 
random direction superposed on the systemic flow ve- 
locity. As the plasma flows across the shock and then 
loses energy as it advects downstream, variations in 
the synchrotron and inverse Compton flux, as well as 
in the linear polarization, occur over a range of time- 
scales. Since the standing shock lies parsecs from the 
central engine, radiation from a hot dust torus and a 
Mach disk at the vertex of the shock cone can both 
be significant sources of seed photons from the scatter- 
ing. Emission from the latter varies along with the en- 
ergy density of the plasma, although the photons from 
an outburst require months to diffuse across the cells; 
this smoothes the fluctuations of the seed photon fleld. 
Traditional synchrotron self-Compton (SSC) scatter- 
ing of the synchrotron radiation emitted by each of 
the cells other than the Mach disk is not included in 
the current, non-parallelized version of the code owing 
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Figure 1: Schematic of the TEMZ model showing the shock, rarefaction, and cylindrical cells. The emission is assumed 
to occur between the conical standing shock and the rarefaction. The only practical restriction on the number of cells 
NcToss across the jet cross-section (as seen in view down axis on bottom) is the memory of the computer, although the 
mean observed degree of linear polarization (Ilopt) at optical wavelengths provides a guideline; Nctobb ^ 1.4(nopt)~^ . 



to the excessive cpu time needed. 

The motivation for the model arises from sahent 
observed properties of the multi-frequency behavior 
of blazars. Most of these properties are ignored as 
"weather" in other models, but are used by the au- 
thor to both shape the TEMZ model and provide 
tests of its abihty to represent the physics of blazars. 
The critical observations include: (1) red- noise power 



l2010llChatteriee et al.l[2008ll2012llLarssonll2012l| .72i 



shorter time-scales of varia bility of flux an d polariza- 
tion at higher frequencies [Marscheiil2012| . (3) mean 
polarization levels as well as fractional deviations from 
the mean that are higher at optical than at lowe r fre- 
quencies (D'Arcang clo 2010:. iJorstad et al.|[2007| . (4) 
apparent rotations in polarization position angle that 
are really just random walks of the projected magnetic 



field direction [P'Arcangelol l201Ct iD'Arcangelo et ahl 
|2007|, (5) breaks in the synchrotron spec trum by more 
than the radiative loss value of 0.5 [e.g.. lWehrle et al.1 
|2012J |. and (6) flares that are often sharply peaked or 
contain multiple peaks, neither of which is reproduced 
by single- or few-zoned models. The dependence of 
items 2-4 on frequency is directly related to the change 
in spec tral index beyond th e break, according to the 



spectra of flux variations in blazars [e.g., Abdo et all model 



Marscher and Jorsta d 2010]. 
An important feature of the model is the differ- 
ent value of the maximum electron energy 70, max 
that is achieved by particle acceleration at the shock 
front. This can either be determined randomly 
from within a po wer-law probability distribution 
Marscher and Jorstad 20101 or as a function of the 



angle that the magnetic field of the cell subtends to 
the shock normal. The simulations presented here 
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adopt 7o,max times cos^(i?||/i?), where the subscript 
II indicates the component of the magnetic field that 
is parallel to the shock normal. The value of 70, max 
is not allowed to fall below some minimum value, 
which is set as an input parameter. Future devel- 
opment can include results from Monte Carlo calcu- 
latio ns of particle acceleration in relativistic shocks 
[e.g., ISummerlin and Barine][201^ . 

The code currently includes synchrotron radiation, 
inverse Compton scattering (IC) of seed photons from 
a hot dust torus, and IC of synchrotron + SSC radia- 
tion from relatively slowly moving plasma in a Mach 
disk at the vertex of the conical shock. The last of 
these is abbreviated as "SSC-MD." If a Mach disk is 
present, this can be the dominant SSC emission, since 
the seed photons are Doppler boosted in the frame 
of the turbulent cells. The combined effects of co- 
spatiality of the emission regions at different wave- 
bands, non-uniform electron energy distribution, dif- 
ferent magnetic field orientations for different turbu- 
lent cells, and light-travel delays, cause correlations 
of variations at pairs of wavebands but with time lags 
and often a lack of one-to-one correspondence of flares. 
This is similar to the observed be havior of blazars 
[e.g.. lMarsche^l2012l IWehrle et all [2012] . The code 
generates simulations of flux (at 68 frequencies from 
lO^'^ to 10^^ Hz) and polarization (at radio to optical 
frequencies) vs. time. 

The synchrotron flux density of a cell depends on 
(among other factors): (1) the volume filling factor 
of electrons with energies high enough to radiate at 
the frequency of observation, (2) the spectral index a, 
determined by the slope of the electron energy distri- 
bution, (3) the normalization factor Nq of the elec- 
tron energy distribution, and (4) the strength B and 
angle ip (corrected for relativistic aberration) of the 
magnetic field relative to the line of sight through the 
factor (Bsinf/')^+". 

Meanwhile, the inverse Compton X-ray or 7-ray 
flux density depends on: (1) the volume filling fac- 
tor of electrons with energies high enough to scatter 
the highest-frequency seed photons to the observed 
photon energy, (2) the slope of the electron energy 
distribution, (3) the normalization factor Nq of the 
electron energy distribution, and (4) the energy den- 
sity of seed photons, which is constant for photons 
emitted by a dust torus and nearly constant if from 
the broad emission-line region, but variable if from a 
relatively slow region inside the jet, such as a Mach 
disk. 

Since some of these factors are different for the two 
processes and for the different sources of seed photons, 
there will be a correlation but not a one-to-one cor- 
respondence between light curves at two widely sep- 
arated frequencies despite the electron energies being 
roughly the same for optical and 7-ray emission. Even 
for relatively closely spaced frequencies (e.g., optical 
and near-IR), differences in volume filling factor as 
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Figure 2: Sample spectral energy distribution produced 
by the TEMZ code, with physical parameters selected to 
produce emission similar to that of BL Lac during a 
7-ray flare (see Fig. [5]). Note the power-law shape in the 
near-IR, optical, and near-UV range with a spectral 
index of 1.3, significantly steeper than both the input 
lower-frequency value of qiow = 0.55 and the 
high-frequency value ctiow +0.5 = 1.05 that one would 
expect solely from radiative energy losses. This results 
from the combination of radiative losses and the 
dependence of the highest energy to which electrons are 
accelerated by the shock on the angle between the 
magnetic field and the shock front (see text). 



a function of maximum electron energy will weaken 
the correlations of the light curves and polarization 
vs. time. At very high photon energies, "orphan 
fiares" are possible if, e.g., the magnetic field of the 
cell that contains the highest energy electrons points 
directly along the line of sight. In this example, no 
synchrotron fiare would be seen even during a major 
7-ray fiare. A Mach disk at the end of the conical rec- 
oUimation shock can provide a variable source of seed 
photons that is well beamed in the frame of the plasma 
passing through the recoUimation shock but not in the 
observer's frame, since the flow speed of the plasma 
beyond the Mach disk is only mildly relativistic. Tem- 
porarily bright regions in a relatively slow sheath of 
the jet could produce a similar effect. The code allows 
each cell to have a randomly oriented turbulent veloc- 
ity relative to the general jet flow. This, along with 
the different maximum electron energy and magnetic 
field direction in the different cells, as well as the small 
volume of each cell, allows for ver y fast variability at 
optic al and 7-ray frequencies [see iNaravan and PiranI 
I2012J . 

If one compares the results of a simulation (Fig. 
51) with actual da ta [see Fig. [S] as well as, e.g. 



Wehrle et al.l 12012 1. the results look quite promising. 
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Figure 3: Sample simulated light curves and linear polarization v s. time generated b y the TEMZ code. Left: Input 
parameters similar to those inferred observationally for 3C 454.3 [Wehrle et a l. 2012], with A'cross ~ 60 (see the caption 
of Fig. [1} . Bottom panel focuses on the brightest fl are seen in the top panel. Right: Input parameters similar to those 
inferred observationally for BL Lac [Marscherll2012l ] . with A'cross = 126. Note the changing ratios of high-energy to 
optical flux, sharply peaked flares, more rapid variations at higher frequencies (7-ray vs. X-ray and optical vs. 
mm- wave), stronger linear polarization at higher frequencies (optical vs. mm- wave), apparent rotations in polarization 
position angle x, and preferred polarization position angle same as jet direction (0° in the simulation). Some examples 
of orphan 7-ray flares are seen. All of these are observed properties of blazars. 



One can see in the plots a number of cases of ap- 
parently "orphan" 7-ray flares, and many flares are 
sharply peaked rather than rounded, although out- 
bursts at lower frequencies are more rounded because 
they occur over larger volumes. Apparent rotations in 
the polarization vector occur in both clockwise (de- 
creasing x) and counterclo ckwise directions, as ob- 
served in some blazars [e.g. iHavashida et al.ll2012| . 

Future development includes parallelization of a 



CH — h version of the code (which was originally writ- 
ten in Fortran-77) so that it can be run efficiently on 
a multi-processor computer. This will allow the inclu- 
sion of SSC of time-delayed seed photons from cells 
other than the Mach disk, thought to be important at 
least in BL Lac objects, if not quasars. In addition, 
emission-line radiation from a cloud that lies nex t to 
the parsec-scale jet [see iLeon-Taveres et al.ll2013J | will 
be added. This can be an additional source of time- 
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51.55 days n= 4.8% x = 



■16.4° 



distinct polarization structure produced by turbulent 
plasma crossing the standing conical shock (see Fig. 
Ill) is seen in nature. 



jet flow 
direction 




Figure 4: Sample polarized intensity map of the emission 
region between the standing shock and the rarefaction 
shown in Fig. [T] at time 51.55 days in the same 
simulation as the one that produced the right-hand side 
of Figure [S] The contours are in factors of 2 starting at 
1% of the maximum in polarization intensity. The 
lengths of the line segments are proportional to the 
polarized intensity, and they are aligned with the 
polarization direction. The orientation is similar to that 
of Fig. [T] with the jet base to the left. The integrated 
degree 11 and position angle x of polarization are 
indicated, with x = 0° parallel to the jet axis (left-right 
in the figure) and positive values of x counterclockwise to 
this direction. Note that, despite the near-cancellation of 
the polarization from different parts of the jet, an image 
of sufficiently high angular resolution would reveal strong 
gradients in both the degree and position angle of 
polarization. 



variable seed photons, since the excitation of the cloud 
will depend on the flux of UV photons emitted by the 
jet. 

The eventual goal is to run the code over a range 
of input parameters to determine the physical condi- 
tions under which various observational features can 
be reproduced. Each simulation is a unique com- 
bination of both random and systematic processes, 
and so will not reproduce actual light curves. How- 
ever, statistical analyses of the timing characteristics 

— correlations across frequencies, power spectra, and 
variations of polarization in both time and frequency 

— can be compared with the same characteristics 
of the data obtained from comprehensive monitor- 
ing programs such as that led by the author's group 

Marscher et al.l [20101 1 . In addition, ultra- high reso- 



lution polarized intensity images, such as those now 
beco ming available with VL BI at a wavelength of 3 
mm [Marti- Vidal et aLll2012 [. can verify whether the 
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Figure 5: Multi-waveband flux and degree P and electric-vector position angle x of optical polarization vs. time for BL 
Lac in 2011. Note rapid variations in P and x, indicative of a turbulent magnetic field. Data are from Marscher et al. 
(in preparation). 
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